Mass Loss Timescale of Star Clusters 



in External Tidal Field 

Ataru Tanikawa, and Toshiyuki Fukushige 

Department of General System Studies, College of Arts and Sciences, 
University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902 
tanikawa@provence. c. u-tokyo. ac.jp 

^ ■ (Received ; accepted ) 

\ Abstract 
Of 

We investigate evolution of star clusters in external tidal field by means of N- 
body simulations. We followed seven sets of cluster models whose central concen- 
tration and strength of the tidal field are different. We found that the mass loss 
timescale due to escape of stars, t m i OS s, and its dependence on the two-body relax- 
ation timescale, t Thji , are determined by the strength of the tidal field. The logarith- 
mic slope [= d\n(t m \ oss ) / dln^^i)} approaches to near unity for the cluster models in 
weaker tidal field. The timescale and the dependence are almost independent of the 
^ | central concentration for clusters in the tidal field of the same strength. In our results, 

! the scaling found by Baumgardt (2001) can be seen only in the cluster models with 
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moderately strong tidal field. 
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1. Introduction 

The escape of stars from clusters is a long-standing problem in the study of dy- 
namical evolution of star clusters. The theory for the escape was first developed by 
Ambartsumian (1938) and Spitzer (1940). It is based on the assumption that distant encoun- 
ters between cluster stars will gradually set up a Maxwellian velocity distribution. Some stars 
have escape velocity of the cluster and they escape. The timescale that the distant encounters 
set up a Maxwellian velocity distribution corresponds to the two-body relaxation timescale 
(Chandrasekhar 1942, Spitzer 1987): 

t r = 0.065 ^ (1) 
nrrrG^ In A 

where n is the number density of stars, m is the mass of the cluster stars, v m is the average 
velocity of the stars, G is the gravitational constant, and A is the Coulomb logarithm. Every 
time the two-body relaxation time elapses, a cluster loses constant fraction of stars which have 
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escape velocity. Therefore, the mass loss timescale of a cluster is proportional to the two-body 
relaxation time. 

This theoretical argument was confirmed by Baumgardt (2001) (hereafter, B01) in the 
absence of the external tidal field. He performed iV-body simulations of evolution of star clus- 
ters, and showed that the mass loss of the cluster happened on the two-body relaxation timescale 
for energy cutoff clusters, and on the two-body relaxation with backscattering correction for 
radial (spatial) cutoff clusters. 

However, the mass loss timescale from the clusters in external tidal field is remained 
unclear. The Collaborative Experiments (Heggie et al. 1998) demonstrated that the mass loss 
timescale does not scale with the two-body relaxation timescale, where multi-mass clusters 
moving in a steady tidal field were simulated. It turns out that the complication is due to the 
existence, at the initial setup, of population of stars that have energies above the escape energy 
and remain in a cluster before finding exits ("potential escapers"). Fukushige, Heggie (2000) 
(hereafter, FH) quantified the escape timescale of the potential escaper as a function of energy, 
and pointed out the escape timescale of the potential escaper is long enough to influence the 
whole mass loss timescale of clusters. 

B01 performed iV-body simulations of star clusters, where equal-mass clusters of W = 
3 King profiles (King 1966) evolve in the steady external tidal field. He investigated the 
dependence of mass loss timescale on two-body relaxation times using simulations whose particle 
numbers are N = 128 to 16384. He found that the mass loss timescale of clusters, t m i OS s; are 
proportional to, t r h 3 ^ 4 , the half-mass relaxation time to the power of 3/4, where the half-mass 
relaxation time is given by 



t th = 0.138 _ * , 2 

vWGln^iV) V ' 



(Spitzer 1987) where N is the number of stars, is the half-mass radius, and 7 = 0.11 (Giersz, 
Heggie 1994). Baumgardt, Makino (2003) found almost the same scaling law for multi-mass 
clusters in time-dependent tidal field. 

However, the t m ioss <x ^rh 3//4 scaling is not convinced as an asymptotic behavior at larger 
N (or t r h) for the following reason. In the limit of large N, the two-body relaxation timescale, 
t r h, is much longer than the escape time delay, t e , which is the duration from the moment when 
the star gets energy above the escape energy to that when it actually escapes from the cluster. 
Therefore, the mass loss timescale should be determined only by the half-mass relaxation time, 

i.e. tmloss ^rh- 

The purpose of this paper is to investigate whether the t m ioss oc ^rh 3//4 scaling obtained by 
B01 is an asymptotic behavior by means of iV-body simulations. We perform sets of simulations 
of star clusters whose two-body relaxation timescale is larger than that of B01 and that are in 
the external tidal field with different strength. 

The structure of the paper is as follows. We describe in detail the model of the cluster in 
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Table 1. Initial cluster models 



Model name 


W 


r t ,i 


n,i/r kg (W = 3) 


r t ,i/r ks (W = 7) 




N 


Edy,r=r t ,i 


k3r0.8 


3 


2.50 


0.80 


0.36 


128- 


- 131072 


7.9 


k3rl.O 


3 


3.13 


1.0 


0.44 


128- 


-131072 


11 


k3rl.3 


3 


4.20 


1.3 


0.60 


128 


- 32768 


17 


k3r2.2 


3 


6.98 


2.2 


1.0 


128 


- 32768 


37 


k3r4.5 


3 


14.0 


4.5 


2.0 


128 


- 16384 


100 


k7rl.O 


7 


3.13 


1.0 


0.44 


128- 


- 131072 


11 


k7r2.2 


7 


6.98 


2.2 


1.0 


128 


- 32768 


37 



section 2. We show the results of iV-body simulations in section 3. section 4 is for discussion. 
In section 5, we summarize this paper. 

2. Simulation Method 

We investigate the evolution of star clusters in the external tidal field by means of TV- 
body simulations. We adopt the conventional model in which the cluster is assumed to move 
on a circular orbit in a spherically symmetric galaxy potential, taken to be that of a distant 
point mass M g . We set the initial center of mass of the cluster at the origin (x,y,z) = (0,0,0), 
with axes oriented so that the position of the galactic center is (— R g ,0,0). We assume that 
the size of the globular cluster is much smaller than R g . If the globular cluster rotates around 
the galactic center at an angular velocity f2 = (0,0, cu), the equation of motion of an individual 
cluster member can be expressed as 

v ■ (Lv ■ 

= -V$ C(i - 2tl x + w 2 (3^e x - Zi e z ), (3) 

where T{ is the position of the i-th particle, and e x , e z are unit vectors that point along the x, 
z axes, respectively. The first term on the right-hand side in equation (3) is the gravitational 
acceleration from other particles in the cluster, the second term is the Coriolis acceleration, 
and the third term is a combination of the centrifugal and tidal forces. We used standard units 
(Heggie, Mathieu 1986), such that M{ — G— —4:E C = 1, where M ; is the initial total mass and E c 
is the initial total energy within the cluster. In our simulations we used a softened gravitational 
potential expressed as 

where rrij is mass of j-th particle and e is a softening parameter. We set the softening parameter, 
e, as 1/32. 

We used King's models (King 1966) to generate initial distribution of star clusters. We 
perform seven sets of simulations of star clusters whose initial dimensionless central potential, 
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Fig. 1. Tidal radii, r tj i , together with the particle distribution (N = 4096), for models with Wo = 3 King 
model. The dashed circles show the tidal radii (r tj i) determined by the external tidal field outward for 
models k3r0.8, k3rl.O, k3rl.3, k3r2.2 and k3r4.5, and the solid circle shows the King cutoff radius (Yk g ) 
beyond which the density is zero in King model. 

Wo, of King model and initial tidal radii determined by the external field, r tj i are different, which 
are summarized in table 1. The number after "k" in model name indicates Wo of King model. 
We set the strength of the external tidal field by giving the angular velocity to in equation (4) 
using a relation, u = ^GMj?>r t ^. The number after "r" in the model name indicates the tidal 
radius scaled by r\ g of Wo = 3, where r^ g is the radius beyond which the density is zero in King 
model (hereafter, King cutoff radius). The values for r^g in the standard unit are 3.13 and 6.98 
for Wo = 3 and 7, respectively. Figure 1 illustrates the tidal radii, r t ,i for simulations of Wo = 3 
King model. The number of particles used for runs are listed in table 1. All particles have 
the same mass, m = M\/N . We perform five runs whose realization of particle distribution are 
different for each TV when TV 5= 8192, and one run which N ^ 16384. 

We perform numerical integrations of equation (3) using a leap-frog integration scheme 
with shared and constant timestep. The stepsize, At, is set to be as 1/64 in models k3r0.8, 
k3rl.O and k7rl.O, and as 1/128 in models k3rl.3, k3r2.2, k3r4.5 and k7r2.2. 

The force calculation is done with the Barnes- Hut algorithm (Barnes, Hut 1986) on 
GRAPE-5 (Kawai et al. 2000), a special-purpose computer designed to accelerate iV-body 
simulations. Actually, we used the same code as in Fukushige, Suto (2001). We use only the 
dipole expansion and the opening parameter 9 = 0.5. It took about 500 CPU hours to complete 
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the most time-consuming run, N = 32768 of model k3r2.2 (about 7.2 x 10 6 timesteps). For 
smaller N simulations, the force calculation is done by direct summations (when N ^ 4096) 
and on host computer (without GRAPE-5, when A" ^ 512). 

Contrary to the standard star cluster simulations, our simulation uses the softened grav- 
itational potential and the leap-flog integrator with relatively large stepsize. And also, the force 
calculation is done with the tree algorithm. These approaches are adopted in order to make the 
two-body relaxation timescale as large as possible. As will be discussed in the appendix, the 
approaches we adopted do not influence the results concerning on the escape from the cluster. 

3. Results 

3.1. Evolution of Total Mass 

Figure 2 shows evolution of the total mass for all cluster models. The curves indicate 
the decrease in mass of the cluster defined by a tidal boundary. We define geometrically the 
cluster member as all stars within the tidal radius from the center of mass of the cluster. The 
tidal radius is expressed as 

-(ST. 

where M is total mass of the members of the cluster at a given time. Since M depends on r t 
itself, some iteration is usually required. We remove stars when they escape far enough (more 
than 4096 from the coordinate origin). 

3.2. Mass Loss Timescale 

Figure 3 shows the mass loss timescale, t m i OS s; of clusters as a function of the initial 
half-mass relaxation time, t^u for Wq = 3 King cluster models. The mass loss timescale is here 
defined as the time when 50 % of the initial total mass is lost. Since we perform five runs when 
N S 8192, we adopt the means of these runs as the mass loss timescale, t m ioss- Table 2 shows 
the maximum deviations among these runs. The maximum deviation is defined as the largest 
difference from the mean. Since we use the potential softening and fix the softening parameters 
for all models, the initial half-mass relaxation time, £ r h,i, is expressed as 

/v>, .3/2 

t Ai = 0.138—^ — — ^ , (6) 

M i 1 / 2 G 1 /21n(0.25r h>i /e) 

where is the initial half-mass radius and e = 1/32. Thus, the initial half-mass relaxation 
timescale, t r h,b is estimated as t^ = 470 x (AT/8192). The initial half- mass radii in the standard 
unit are r h i = 0.84. 

As can be seen from figure 3, a cluster in stronger tidal field lose its mass sooner. The 
difference of the mass loss timescale is considered to be due to the difference of the escape 
energy (the potential height at the Lagrangian point), which is lower for a cluster in stronger 
tidal field. Figure 4 shows change in energy of an individual star for k3rl.O and k3r2.2 models 
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Fig. 2. Evolutions of the mass of the clusters. The numbers under the curves indicate the numbers of 
particles, N. 



Table 2. Maximum deviation of t m i oss for N ^ 8192. 



N 


k3r0.8 


k3rl.O 


k3rl.3 


k3r2.2 


k3r4.5 


k7rl.O 


k7r2.2 


128 


11 


32 


62 


40 


86 


41 


151 


256 


21 


20 


41 


65 


145 


33 


94 


512 


22 


19 


40 


90 


129 


48 


53 


1024 


50 


108 


40 


38 


189 


57 


121 


2048 


37 


40 


115 


230 


138 


36 


123 


4096 


76 


81 


145 


194 


419 


51 


183 


8192 


92 


166 


201 


118 


286 


200 


211 
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Fig. 3. Mass loss timescale of the clusters as a function of initial half-mass two-body relaxation timescale 
for Wo = 3 King cluster models. The numbers near the curves indicate rt,i/rkg- 

(N = 1024). In figure 4 we plot {-Emax,i}med, the median value of the maximum energy records 
at a given time, t\] 

£max,i(*i) = max-j^t)}, (7) 

t<ti 

1 3 1 

Ei{t) = ~{vj + v y 2 + vj) + $ c ,i + (--x 2 + -z 2 ) (8) 

of randomly selected 50 particles, where v x i, v y i, and v z i are x, y, and z components of velocity 
of i-th particle, respectively. We do not update E mauXti (ti) after stars escape from the cluster. 
The representative value, {-Emax,i}med, may trace the energy acquired by two-body relaxations, 
on average. When stars obtain enough energy, they escape from cluster eventually. The dashed 
lines indicate the escape energy of each models, expressed as 
3 GM 

£ cri t = -- — • (9) 

As shown in figure 4, the increase of {-Emax^jmcd at t < 600 is very similar among models 
k3rl.O and k3r2.2, which indicates that the relaxation processes in both models occur with 
almost the same timescale. However, because the escape energy in weaker tidal field is larger, 
it takes more time for a significant fraction of stars to escape in model k7r2.2 than in model 
k7rl.O. 

Figure 5 shows logarithmic slopes, a(t r h,i) = (iln(t mloss )/(iln(t r i lji ), of the dependence on 
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Fig. 4. Change in energy of an individual star for k3rl.O and k3r2.2 models (N = 1024). The dashed and 
solid lines indicate {-E m ax,i}mcd (the definition is in the text) in models k3rl.O and k3r2.2, respectively. 
The upper dotted line indicates the escape energy in model k3r2.2 and the lower dotted line indicates that 
in model k3rl.O. 

t r h,i of the mass loss timescale shown in figure 3. Actually, we calculated the logarithmic slope 
«(^rh,i) using a relation; 



where t m ioss(^rh,i) is the mass loss timescale when the initial half-mass relaxation timescale is 



As can be seen in figure 5, the logarithmic slope depends on the the strength of external 
tidal field. In model k3rl.O, where the initial distribution (Wo = 3 King model) and the tidal field 
( r t,i = r k g ) are set to be the same as in B01, we reproduce the asymptotic power, a(t r h,i) ~ 0.75, 
obtained by him. Note that the initial half-mass relaxation timescale of the run with the largest 
N in B01 is about 230. On the other hands, when the tidal field is stronger (model k3r0.8) the 
slope is smaller [a(t r h,i) ~ 0.4], and when it is weaker (models k3rl.3, k3r2.2, and k3r4.5) the 
slope is larger than a(t r h,i) = 0.75. The slope a(t r h,i) for the models in much weaker tidal field 
(more k3r2.2 and k3r4.5) seems to approach asymptotically a(t r h,i) ~ 0.9, not 1. 

The dependence of the slope is considered to be associated with the population of the 
potential escaper. For the clusters whose initial half-mass relaxation timescale is smaller (or 
smaller N), the escape time delay due to the potential escaper delays more the total mass loss 



a(£rh,i) 




(10) 
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Fig. 5. Logarithmic slope a(t r h,i) [ = <iln(i m ioss)/rfln(Ah,i)] of dependence on the initial half-mass relax- 
ation timescale. The dashed line indicates a(t r h,i) = 0.75. 

of clusters. Figure 6 shows the evolution of the fraction of potential escapers M pe /M in the 
King Wo = 3 cluster models (N = 16384, except for model k3r4.5, N = 8192). We can see that 
as the slope a(t r h,i) shown in Figure 5 decreases, the fraction M pe /M increases. 

3. 3. Effect of Initial Concentration 

We investigate the dependence of the mass loss timescale on initial concentration of the 
cluster profile with W = 3 and 7 King profiles. Figure 7 shows the mass loss timescale, t m i OS s, 
as a function of the initial half-mass relaxation time, t r h,i, for models k3rl.O, k3r2.2, k7rl.O and 
k7r2.2. In the models k7rl.O and k7r2.2, the initial half-mass relaxation timescale is estimated 
as t r h,i = 420 x (iV/8192). The initial half-mass radii in the standard unit are r^i = 0.77, for 
Wq = 7. As seen in figure 7, the differences between models k3rl.O and k7rl.O and between 
models k3r2.2 and k7r2.2 are small, which means the mass loss timescale does not depend 
so much on the concentration of the cluster. Figure 8 shows the logarithmic slope, a(t r h,i) [ 
= d\n(t mloss ) / dln(t rhti ) ). 
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Fig. 6. Evolution of the fraction of the mass of potential escapers in the King Wo = 3 clusters. 
4. Discussion 

4-1. Mass Loss Timescale for Small N 

We discuss on the small t r h,i (or small N) limit of the mass loss timescale, t m i OS s, of 
clusters. When N is small, the half-mass relaxation time, t r h,i, is sufficiently small compared 
with the escape time delay, t e , and it is expected to dominate the total mass loss timescale, 
^mioss- In figure 9, we show the mass loss timescale of clusters scaled by the average dynamical 
time within the tidal radius, t dyjr=rti , which is expressed as 

W= rtli = G1/ 2 M i/2 C 11 ) 

and listed in table 1, as a function of the initial half- mass relaxation timescale for Wo = 3 King 
cluster models. The escape time delay is proportional to the average dynamical time (FH). In 
this figure, we don't plot that for model k3r4.5, since only this model with N = 128 and 256 
experiences core collapse before the half mass is lost. As can be seen in figure 9, all of the 
scaled mass loss timescales, t r h,iAdy,r=r t v approach to ~ 10 around t r h,i ~ 7, which means that 
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Fig. 7. Mass loss timescale as a function of the initial half-mass relaxation timescale for models k3rl.O, 
k3r2.2, k7rl.O and k7r2.2. 




Fig. 8. Logarithmic slope a(t r h,i)[= dln.(t m \ oss ) / d]n(t T h. t i)) as a function of the initial half-mass relaxation 
timescale for models k3rl.O, k3r2.2, k7rl.O and k7r2.2. The dashed line indicates a(t r h,i) = 0.75. 
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Fig. 9. Mass loss timescale scaled by the averaged dynamical time within the tidal radius as a function 
of the initial half-mass relaxation times for Wo = 3 King cluster models. 

the mass loss timescale is determined mainly by the averaged dynamical timescale or escape 
time delay in small t r h,i (or small N) limit. 

4-2. The Asymptotic Slope: a(t r h,i) ~ 0.9 

In section 1, we conjecture that the mass loss timescale should be determined only by 
the two-body relaxation timescale at its large limit. Our simulation results show that when 
the external tidal field is not strong (ex. model k3r2.2) the mass loss timescale is almost 
proportional to the half-mass relaxation timescale. However, its logarithmic slope seems to 
converge to slightly small value, a(t r h,i) ~ 0.9. Here, we discuss origins of the small value 0.9. 

One possible answer is that the two-body relaxation timescale we used may not correctly 
scale the real relaxation process occurred in the cluster. Another possible answer is that our 
simulation time span may be still not long enough to exclude the influence of the escape time 
delay on the mass loss timescale. Our simulation result and analytical estimate appear to favor 
the later possibility. 

Figure 10 shows evolution of the Lagrangian radii as a function of time scaled by the 
initial half-mass relaxation time, t r h,i for model k3r2.2 with N = 1024,2048,4096. We can 
see that the 20% Lagrangian radii of clusters with different N are in good agreement, which 
indicates that the half-mass relaxation time is surely a good measure of the internal relaxation 
process. 
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Fig. 10. Evolution of the Lagrangian radii containing 20, 50 and 70 % of the total mass as a function of 
time scaled by the initial half-mass relaxation time. 

We estimate distribution of the escape time delay in model k7r2.2 using the results of 
FH, which gave the distribution of escape time delay for fixed King potential in r tii = 7\ g tidal 
field. According to this paper, the fraction of the initial potential escapers staying in a cluster, 
P(t), is expressed by 

P(t) =&^{[l - Uon(E)} (1 + 3.96^£ V' 729 + f non (E)\, (12) 

E 

where E — \{E — E crit ) / E CTit | and E is the energy of a potential escaper, g dis (E) is the fraction of 
the initial potential escapers, and f n0 n{E) is the fraction of non-escaper, which can not escape 
from the cluster due to the regular orbits, in the initial potential escapers. We took 0.04, 0.06, 
0.08, 0.12, 0.16, and 0.24 as representative values of E. 

Figure 11 shows the distribution of escape time delay equation (12). In this figure we 
can see that at t ~ 10 4 , which corresponds to the maximum of our simulation span, only half 
of the initial potential escapers have been depleted. This means that the mass loss timescale, 
^mioss, of clusters are still affected not only by the initial half-mass relaxation time, t r h,i, but by 
the escape time delay, t c within our simulation span. At least, 10 times simulation span may 
require to exclude the influence of the escape time delay. 
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Fig. 11. Evolution of the fraction of the initial potential escapers based on the results of FH in the case 
of k7r2.2 series. Dashed line indicates t = 40000, which is equal to the largest simulation span. 

5. Summary 

We investigate the evolution of star clusters in external tidal field by means of N- 
body simulations. We follows evolution of seven sets of clusters, whose dimensionless central 
potential, Wo, of clusters and the strength of external tidal field are different. 

Our main conclusions are the following: 

1. The mass loss timescale depends on the strength of external tidal field. 

2. The Logarithmic slopes, a(t rh i ) [ = <iln(t mloss )/<iln(t rhji )], also depend on the strength of 
external tidal field. In model k3rl.O, whose run parameters are set to be the same as in 
B01, we can reproduce the asymptotic power, a(i r h,i) ~ 0.75. However, when the tidal 
field is stronger (model k3r0.8) the slope is smaller [a(t r h,i) ~ 0.4], and when it is weaker 
(models k3rl.3, k3r2.2, and k3r4.5) the slope is larger than a(t r h,i) =0.75. The slope a(£ r h,i) 
for the models in much weaker tidal field (more k3r2.2 and k3r4.5) is seen to approach 
asymptotically to near unity, but exactly a(t r h,i) ~ 0.9. 

3. The mass loss timescale is almost independent of the dimensionless central potential, Wo- 

We are grateful to Jun Makino and Holger Baumgardt for many helpful discussions. This 
research was partially supported by the Grants-in-Aid by the Japan Society for the Promotion 
of Science (14740127) and by the Ministry of Education, Science, Sports and Culture of Japan 
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Fig. 12. Evolution of the total mass of cluster in model k3r2.2 (N = 1024). and solid lines are for runs 
with and without GRAPE-5, respectively. 

(16684002). 

Appendix 1. Reliability of Simulations 

A. 1.1. Force Accuracy 

We use GRAPE-5 and the tree algorithm for the force calculations. Both arise artificial 
errors in the calculated force. When calculated with GRAPE-5 the force between two particles 
has a relative error of ~ 0.2 percent, because of the low accuracy of internal expressions in 
the GRAPE-5 chip (Makino, Ito, Ebisuzaki 1990, Kawai et al. 2000). The force calculation by 
the tree algorithm also contains error of similar magnitude, though it depends on the opening 
parameter. However, such errors caused by GRAPE-5 and the tree algorithm may not influence 
so much our results, because they are small compared to two-body relaxation effects (Hernquist, 
Makino, Hut 1993), which we want to handle correctly in our simulations. Figure 12 shows 
evolution of the total mass of clusters when force calculations are done by the direct summation 
with and without GRAPE-5, model k3r2.2 (N = 1024). The difference between runs with and 
without GRAPE-5 is smaller than its run-to-run variation. Figure 13 shows evolution of the 
total mass calculated using direct summation and the tree algorithm (6 = 0.5) in model k3r2.2 
(N = 8192). Both runs used GRAPE-5. The difference of the results is sufficiently smaller than 
its run-to-run variation. 



15 



i 1 1 1 1 1 r 




j i i I i i L 



8000 16000 

t 

Fig. 13. Evolution of the total mass of clusters in model k3r2.2 (N = 8192). The solid and dashed lines 
indicates those obtained using direct summation and the tree algorithm, respectively. Three runs whose 
initial realization are different are performed for each case. 

A. 1.2. Accuracy of Time Integration 

In our simulations, we use the leap-frog scheme, despite of the dependence of forces 
on the velocities. It is well-known that the nature of leap-frog scheme becomes worse in the 
presence of this dependence. However, figure 14 shows evolutions of the total mass of clusters 
when At = 1/128 and 1/256 are converged, which means the leap-frog scheme with At = 1/128 
does not influence our results. 

A. 1.3. Potential Softening 

We use the potential softening for the force calculation. One might suspect if the two- 
body relaxation doesn't happen in such smoothed potential. Using the potential softening, 
close encounter between particles is surely suppressed. However, energy of particles is changed 
dominantly by accumulation of distant encounters, not by a few close encounters (Spitzer 1987). 
So, two-body relaxation should occur in the smoothed potential. 

We check that the two-body relation happens with the softened potential. The upper 
panel of figure 15 shows the evolutions of the minimum potential in the clusters, $ m i n , as a 
function of time scaled by the initial half-mass relaxation time, t r h,i, in the clusters of model 
k3r2.2 of N = 1024,2048,4096,8192 and e = 1/32. The central potential decreases owing to the 
two-body relaxation effect. As seen in the upper panel of figure 15, these curves are in good 
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Fig. 14. Evolution of the total mass of clusters in model k3r2.2(iV = 8192). The solid and dashed lines 
indicates those obtained using timcstcp size At = 1/128 and At = 1/256, respectively. Three runs whose 
initial realization are different are performed for each case. 

agreement and evolution of the central potential is scaled by the two-body relaxation timescale. 

The lower panel of figure 15 shows the evolutions of $ m i n in the clusters of model k3r2.2 
of N = 8192 and e = 1/32,1/48,1/64 as a function of timescale by the two-body relaxation 
timescale [equation (6)]. These curves are also in good agreement, which indicates that even 
with a softened potential, two-body relaxation occurs on the half-mass relaxation timescale, 
described by equation (6). 
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Fig. 15. Evolution of the minimum potential in clusters as a function of time scaled by the initial half-mass 
relaxation time in model k3r2.2 of N = 1024,2048,4096,8192 and e = 1/32 (top), and N = 8192 and 
e = 1/32,1/48,1/64 (bottom). 
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